Hydraulic tortuosity in arbitrary porous media flow 
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Tortuosity (T) is a parameter describing an average elongation of fluid streamlines in a porous 
medium as compared to free flow. In this paper several methods of calculating this quantity from 
lengths of individual streamlines are compared and their weak and strong features are discussed. An 
alternative method is proposed, which enables one to calculate T directly from the fluid velocity field, 
without the need of determining streamlines, which greatly simplifies determination of tortuosity 
in complex geometries, including those found in experiments or 3D computer models. Numerical 
results obtained with this method suggest that (a) the hydraulic tortuosity of an isotropic fibrous 
medium takes on the form T = 1 +py/l — ip, where ip is the porosity and p is a constant and (b) the 
exponent controlling the divergence of T with the system size at percolation threshold is related to 
an exponent describing the scaling of the most probable traveling length at bond percolation. 



I. INTRODUCTION 

Investigation of transport through porous media is of 
paramount importance in many areas of science and en- 
gineering. One of the main problems is to find out how 
the value of permeability, which synthetically describes 
how the flow is retarded by the porous medium struc- 
ture, can be related to some more fundamental, well- 
defined parameters determined solely by the geometry of 
the medium, as such relation could be used, for example, 
to fabricate materials of desired physical properties. 

One of the most well-known theories of this kind was 
developed by Kozeny and later modified by Carman [![. 
In their approach a porous medium is assumed to be 
equivalent to a bundle of capillaries of equal length and 
constant cross-section. These assumptions lead to the 
semi-empirical Kozeny-Carman formula [3-0] 
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which relates the permeability (k) to four structural pa- 
rameters: the porosity <p>, the specific surface area S, the 
hydraulic tortuosity T, and the shape factor /?. In this 
equation ip is a dimensionless quantity defined as the frac- 
tion of the porous sample that is occupied by pore space 
(0 < ip < 1), S equals to the ratio of the total interstitial 
surface area to the bulk volume, ft is a constant charac- 
teristic for a particular type of granular material, and T 
is a dimensionless parameter defined as 



T : 



(A) 
L 



> 1, 



(2) 



where (A) is the mean length of fluid particle paths and L 
is the straight-line distance through the medium in the 
direction of macroscopic flow. Eq. ([1]) has been found 
to agree well with experimental results for random packs 
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FIG. 1. (Color) Normalized velocity field (u/u m ax) for two 
highly porous systems, ip = 0.95 (left) and ip — 0.99 (right). 
Lighter (darker) colors represent larger (smaller) fluid veloci- 
ties. 



of monodisperse granules, e.g. spheres or well sorted and 
rounded sands. 

The Kozeny-Carman approximation of a porous 
medium can be used to model also other types of trans- 
port, e.g. diffusion or electric current. This observation 
resulted in introducing several distinctive, experimen- 
tally measurable quantities that in the capillary approx- 
imation can be readily linked with the tortuosity of the 
capillaries as given by Eq. @. In this way the term 
'tortuosity' has been overloaded with several essentially 
different meanings [|| , as depending on the context it can 
refer not only to hydraulic, but also geometric @, @|, dif- 
fusional 043, electrical [ 1 , fToMT^ ] . thermal [l3j], acoustic 
[l2| , or streamline [TJ, [l5[ tortuosity, with no clear con- 
sensus on the relation between them. Besides, in the 
literature different quantities, including T -1 , T~ 2 , and 
T 2 have also been called 'tortuosity'. 

While the permeability in the Kozeny-Carman theory 
depends on four structural factors: porosity, specific sur- 
face area, a shape factor and a hydraulic tortuosity, until 
recently only the first two of them could be measured in 
non-trivial cases, and only porosity could be measured 
relatively easily. Carman himself attempted to estimate 
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hydraulic tortuosity by injecting dye into a bed of glass 
spheres and concluded that T ss \[2 [1]. However, in 
practice only the product f3T 2 (known as the Kozeny 
constant) could be determined in non-trivial experimen- 
tal setups, which left the shape factor and hydraulic tor- 
tuosity as essentially indeterminate quantities. This ob- 
servation led several researchers to ponder whether hy- 
draulic tortuosity really exists as a fundamental attribute 
of the pore space or whether it is just a 'fudge factor', an 
adjustable parameter used to fit the model to the exper- 
imental data [1, [HI. In spite of these difficulties, diffu- 
sional and electrical tortuosities are one of the basic pa- 
rameters commonly used to characterize real porous me- 
dia in such diverse areas as medicine [TH, [l?} , marine biol- 
ogy or advanced materials [l8[ , and the hydraulic tor- 
tuosity remains a key concept of many advanced theories, 
e.g. the Effective Medium Approximation 0, [l9j]. More- 
over, modern technology has made it possible to deter- 
mine the fluid velocity field in quite complex geometries 
both experimentally jMlH and numerically [jij Silli- 
ly . This, at least in principle, enables one to determine 
flow streamlines and hence renders the hydraulic tortu- 
osity a measurable quantity. 

At this point, however, there appears an unexpected 
problem. A textbook recipe requires to calculate the hy- 
draulic tortuosity as an 'average path of fluid particles 
through the porous medium' 0] without specifying what 
sort of averaging is actually meant. For example, is it to 
be taken over the whole volume or over a cross-section, 
and in the latter case — are all cross-sections equivalent? 
Should the average be weighted and how? Ambiguity in 
the definition of T was noticed already by Bear Q , who 
remarked that the average pathlines could be obtained 
either by averaging the pathlines of all fluid particles 
passing a given cross-section of the medium at a certain 
instant of time (geometrical approach) or during a given 
period of time (kinematical approach). Bear himself pre- 
ferred the geometrical approach, but he never explained 
how his tortuosity tensor, which is a macroscopic quan- 
tity, could be calculated from microscopic flow stream- 
lines. Clennell [4j gave several convincing arguments in 
favor of the opinion that the hydraulic tortuosity should 
be calculated as a kinematical average in which the path- 
lines are weighted with fluid fluxes. However, until very 
recently a lack of precision in the definition of T was not 
regarded as a problem, as this quantity was considered 
to be too difficult to be calculated in a general case, and 
most attempts in this direction concentrated on rather 
simple models where the results did not depend on the 
averaging procedure. Under these circumstances, when 
in recent years it became possible to simulate flows nu- 
merically with unprecedented accuracy in complex ge- 
ometries in which fluid fluxes continuously change in sec- 
tional area, shape and orientation as well as branch and 
rejoin, researchers developed their own methods of calcu- 
lating T flil. [TEl. l24l - [2l| . Closer inspection of these papers 
leads to a rather surprising conclusion that no consensus 
has been reached as to the actual meaning of the numer- 



ator in Eq. ([2]), as each research group interpreted the 
average in Eq. ^ in their own, unique way. 

The aim of this study is to propose a universal, efficient 
method of calculating hydraulic tortuosity in an arbitrary 
geometry. To this end in Sec. [H] we analyze several algo- 
rithms used so far and show their weak and strong fea- 
tures. Then, in Sec. Mil we present a new method, which 
enables one to calculate T without the need of deter- 
mining any streamlines, which is often an ill-conditioned 
numerical problem especially if only approximate 
values of the velocity field are available. This is the main 
result of the paper. Its significance lies in that it greatly 
simplifies determination of hydraulic tortuosity in exper- 
iments or three-dimensional numerical simulations. In 
Sec. II VI we apply this method in two cases that cannot be 
efficiently treated with other methods: very high porosity 
(fibrous medium) or very low porosity (system at perco- 
lation). Finally, Sec.|V]is devoted to conclusions. 

II. COMPARISON OF EXISTING METHODS 

As it was already mentioned, several methods of calcu- 
lating the hydraulic tortuosity in an arbitrary geometry 
are available in the literature. Most of them reduce the 
problem to calculating T as a weighted average of the 
form 

T (3) 

where i enumerates discrete streamlines, Aj is the length 
of the i-th streamline, and Wi is a weight. 

Knackstedt and Zhang [2(| |27| used Eq. (J3j) with i 
running through the nodes of a regular lattice on the 
inlet cross-section and chose the weights of the form 

1 

w l = -, (4) 

H 

where ij is the time in which a fluid particle moves along 
the i-th streamline. The intention behind (jlj was to 
weight each streamline length proportionally to the over- 
all volumetric flow associated with this streamline. Since 
the streamlines sample the inlet plane uniformly, the 
weights satisfying this condition should be proportional 
to the components of the fluid velocities parallel 

to the macroscopic fluid flow direction (here assumed to 
be directed along the x axis) and measured at the points 
where the i-th streamline cuts the inlet plane. Note that 
the weights defined by ^ might be equivalently written 
as Wi — Vi/Xi, with vi being the average fluid velocity 
along the i-th streamline. However, apart from some triv- 
ial geometries (e.g. straight capillaries of equal length), 
fluid velocity along a typical streamline in a complex ge- 
ometry, for example in a granular porous medium, can 
vary by several orders of magnitude [f4j and hence there 
is no connection between Uj/Aj and (vi)™. For this reason 
it is not clear what the quantity calculated by Knackstedt 
and Zhang has to do with the actual hydraulic tortuosity. 
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Koponen et al. [24| introduced two families of tortu- 
osities T?, and T^f, n £ Z, defined through 



(^) 



and 



(T, 



J A X n (r)v(r) <P 
J A v(r) d 2 r 



J v X n (r)v(r) d 3 r 
J v w(r)d 3 r 



(5) 



(6) 



where A is an arbitrary cross-section perpendicular to the 
macroscopic fluid flow direction, V is the volume of the 
sample, A(r) = X(r)/L is the tortuosity of the flow line 
passing through a point r, v(r) = |v(r)| is the fluid speed 
at r, and v(r) = inside the solid phase. The index 'S' 
at indicates that this quantity is to be calculated on 
a surface (cross-section), whereas 'V at T% indicates a 
volumetric quantity. 

Using a simple capillary model Koponen et al. con- 
cluded that = but it is not difficult to show that 
this is not true in a general case. To this end it suffices to 
consider a bundle of straight and wavy cylindrical cap- 
illaries of the same radius and different lengths. The 
contribution of each of such capillaries to the integrals in 
([5]) is proportional to the area of its cross-section with A, 
which readily implies that is A-dependent. For this 
reason should not be used to calculate the hydraulic 
tortuosity unless, perhaps, in very large, homogeneous 
systems where the effects of A-dependence could be av- 
eraged out. 

In their actual calculations Koponen et al. used Eq. © 
with the integrals approximated by sums [24j . 
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where are some points that sample uniformly the avail- 
able pore space, either by being chosen at random [24| 
or by being identified with the nodes of a lattice used 
to model the system (25|. While for a given steady ve- 
locity field T^y is a mathematically well-defined quantity 
to which the r.h.s. of (0 should converge as the num- 
ber of points Yi goes to infinity, the fact that it is de- 
fined through volumetric integrals introduces some ad- 
ditional, presumably unintentional weighting that favors 
longer streamlines. Using a method described in Sec. IIIII 
it can be shown that for flows of incompressible fluids 
without eddies, T% can be expressed in terms of surface 
integrals 



f A \(r)v ± (r)d*r 



(8) 



where v± is the component of the fluid velocity normal 
to surface A. Comparison of this formula with Eq. ([5]) 
confirms that ^T$. 



Another approach was proposed by Matyka et al. [Tj| , 
whose formula for the hydraulic tortuosity (here denoted 
by T M ) in an integral representation can be written as 



t m = l4 A ( r H( r K r 



J A v ± (r)d 2 r 



(9) 



where the surface A need not be perpendicular to the 
macroscopic flow direction and can even be curved, and 
A and v± are assumed to vanish inside the solid phase of 
the medium. In contrast to Eqs. ([5]) and ([6]), in which 
tortuosities of individual streamlines were weighted with 
local fluid speeds, in Eq. © they are weighted with lo- 
cal fluxes. This guarantees that for incompressible flows 
both integrals in Eq. ©, and hence T M , are independent 
of A. The actual numerical calculations were performed 
using a two-dimensional model system and Eq. ^ was 
approximated with an arithmetic mean 
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where are some points on a cross-section satisfy- 
ing a constant-flux constraint between two neighboring 
streamlines and N is the number of these points. 

Of the three methods presented in this overview, only 
that of Matyka et al. correctly addresses the problem of 
recirculation zones — their contribution to T M vanishes. 
Another advantage of Eq. (fTU)) is that all terms in the 
sum are of the same order of magnitude, whereas the 
sums in ([7]) contain terms that can differ by several or- 
ders of magnitude. It is a consequence of a fact that 
Eq. ([7]) implicitly divides the space into regions of ap- 
proximately equal volume and assigns to them equal im- 
portance, whereas flow in a porous medium takes place 
mainly in a few conducting channels which occupy only 
a small fraction of the porous space. For this reason 
one can expect that for the same number of streamlines 
Eq. (|10|) . which assigns equal importance to equal fluid 
fluxes, will be loaded with a much smaller numerical er- 
ror. A disadvantage of Eq. (|T0|) is that it would be diffi- 
cult to apply it to three-dimensional problems, the main 
difficulty being to find the points satisfying a constant- 
flux condition. 

Note also that it follows from Eqs. ([8]) and ([9]) that 



T 



M 



fll) 



for arbitrary incompressible flows. 



III. ALTERNATIVE METHOD 

Just as it is possible to express a volumetric integral 
© as a surface integral ((5} , it is possible to express a sur- 
face integral Q as a volumetric integral. The resulting 
formula reads 



jiM _ 



J v v x (r) tfV 



(12) 
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where v x denotes the velocity component parallel to the 
macroscopic flow direction. This equation can be written 
in a particularly simple form 



(13) 



in which (...) denotes a spatial average over the pore 
space. This is the main result of the paper, but before 
we proceed to prove it, a few remarks are in order. 

Fluid streamlines are determined numerically by 
solving — usually thousands of times — an ordinary differ- 
ential equation of motion of a massless fluid particle in 
a given velocity field. This field is obtained by some 
method of the computational fluid dynamics (CFD), 
which yields approximate values of the fluid velocity field 
at limited number of discrete points (a mesh). To com- 
pute a streamline, the velocity at an arbitrary point is 
thus required, and this is calculated using some kind of 
extrapolation from the values at the mesh, which in- 
troduces additional errors. For these reasons numeri- 
cal determination of streamlines is time-consuming and 
error-prone. This task is especially difficult at low porosi- 
ties, as in this case the conducting channels are narrow 
and particularly sinuous, and under such conditions the 
computer-generated streamlines tend to hit the channel 
walls and terminate [lj, |24| . The main advantage of (fT2"|) 
is that it enables to calculate a tortuosity without any 
need to determine individual streamlines. It can be thus 
used in very complex geometries, e.g. close to the perco- 
lation threshold in 2D flows or in large-scale, low-porosity 
3D flows. 

The idea that the ratio (v)/(v x ) can be related to the 
hydraulic tortuosity is not new. Carman [l[ used it to 
argue that permeability must be proportional to T~ 2 , 
cf. Eq. ([T]), rather than to T _1 , as had been earlier pos- 
tulated by Karman. Koponen et al. [24j used this ra- 
tio explicitly as one of several possible definitions of the 
hydraulic tortuosity. However, all these attempts were 
based on a simple model where a porous medium is as- 
sumed to be equivalent to a group of parallel channels 
and no attempt was made towards justification of this 
approach in more general cases. 



A. Proof of equation (|12[) 

The system volume V can be divided into two disjoint 
subsets, the porous space V and the solid phase Vb (see 
Fig. [2]). Similarly, any cross-section A can be divided into 
A' = A n V and A = A n V Q . We assume that v = 
and A = at any r £ Vq. Let V* C V denote the set of 
all points r £ V such that the streamline cutting r joins 
the inlet and outlet surfaces. A flow for which V* ^ V 
shall be called reentrant. 

Assume that the flow is stationary, incompressible and 
not reentrant. Incompressibility of flow implies that for 
any cross-section A perpendicular to the flow direction 




FIG. 2. Quantities used in the proof of equation (|12|) . Solid 
phase (Vb) is immersed in porous phase (V). A cross-section 
A consists of the solid (Ao) and porous (A') part. V* is made 
up by all streamlines connecting the inlet and outlet planes. 
Formation of eddies, e.g. in cavities, would violate Eq. (fl2|) . 



the denominator in Eq. © is equal to the total flux 
through the porous sample. This leads to 

/ v x {r)d 3 r = L [ v ± (r)d 2 r. (14) 

JV J A 

To prove (fT2")) . it thus suffices to show that 

/ v(r)d 3 r= [ A(r)t)_L(r)d 2 r. (15) 

JV J A' 

Since v(r) is defined and continuous at each r, and 
A(r) is defined and continuous at each r except for some 
points from a zero-measure subset D C V [14], both 
integrals in (IT51) exist. 

Let A be a cross-section perpendicular to the flow di- 
rection (i.e. to the x-axis) such that each streamline cuts 
A 1 only once (e.g. A is the inlet plane). Let s(r) be 
the distance from A' to r along the streamline passing 
through r. Except for a zero-measure set, any point in V 
can be uniquely identify by the streamline it belongs to 
and s(r). Each streamline, on the other hand, is uniquely 
identified by r £ A' belonging to this streamline. Thus 
A' and s can be used to change the integration variables 
in the integral on the left-hand-side (l.h.s.) of Eq (fTo) 
from x,y,z to s,y,z. Since dx/ds — v x (r)/v(r), and 
J J v x dydz is constant along streamlines in incompress- 
ible flows, we arrive at 



v(r) d 3 r 



v x (y, z, s) dydzds 
ds(r) I v x (r)d 2 r 



\(r)v x (v)d 2 r, 



(16) 



which is the right-hand-side (r.h.s.) of (fT5"]) . Validity of 
(fTS"]) can be extended to cross-sections of arbitrary shape 
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by noticing that v±(r)d 2 r is the flux associated with a 
streamline that passes through r and hence is constant 
along a streamline since the fluid is incompressible. For 
the same reason, if a streamline cuts A many times, con- 
tributions to the r.h.s. of (IT5|) from each subsequent cut- 
ting are the same in magnitude but of alternate signs and 
cancel out in pairs, so that effectively each streamline 
contributes to the integral as if it cut the cross-section 
only once. 

It is interesting to notice that using the same argu- 
ments Eq. (IT51) can be generalized to 

f(r)v(r)d 3 r= f /(r)A(r)u x (r) d 2 r, (17) 

V J A' 

where a function /(r) has a constant value along each 
streamline. 



B. Conditions of applicability of equation (112[) 



Validity of Eq. (fT3| is based on two assumptions: the 
fluid is incompressible and the flow is not reentrant. The 
latter condition is met, for example, for irrotational or 
potential flows. This implies that Eq. (fTB"]) can be used 
to calculate a hydraulic tortuosity for inviscid fluids. An- 
other important class of potential flows are those gov- 
erned by the Laplace equation, e.g. tracer diffusion or 
electric current. Thus we conclude that Eq. (TT3|) can be 
used to calculate diffusional or electrical tortuosities. 

Real fluids, however, are viscous and as they flow 
through a porous medium, some recirculation zones (ed- 
dies) are produced by rapid changes in pore aperture or 
blind pore spaces. These eddies make the flow reentrant 
even at very low Reynolds number. The contribution 
from reentrant zones to the volumetric integral in (1151) is 
strictly positive, whereas it vanishes for the surface in- 
tegral. Therefore, from a mathematical point of view, a 
weaker relation replaces (|15[) for general laminar viscous 
incompressible flows, 



rpM 



(18) 



However, in flows through porous porous media at low 
Reynolds number the volumes where the flow is reen- 
trant not only constitute a small fraction of the total 
porous volume V' , but the fluid velocity in these vol- 
umes is at least an order of magnitude smaller than that 
in conducting channels. Therefore the contribution from 
reentrant regions to the volumetric integral in (|15|) can 
be expected to be negligible and hence Eq. (jT5)l should 
be also applicable to viscous flows, at least in the low 
Reynolds number regime. 



IV. APPLICATIONS 

To verify usability of Eq. (| L3[) , we employed it to find 
the hydraulic tortuosity in a model of freely overlapping 
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FIG. 3. (Color) Tortuosity (T) as a function of porosity (ip) in 
the model of overlapping squares for a = 10 l.u. and L = 1000 
l.u. Symbols represent the results obtained with Eq. (|12|) . 
The vertical dashed line shows the percolation threshold ip c ~ 
0.367 below which no flow (in the limit of L — > oo) is possible. 
Inset: a log-log plot of T — 1 as a function of 1 — ip; the solid 
line is the best-fit to T — 1 oc y/1 — <p for <p > 0.8. 



squares [3, [IH [24|, [HI- In this model one considers a 
two-dimensional lattice with a porous matrix modeled 
with freely overlapping solid squares of size a x a lat- 
tice units (l.u.) placed uniformly at random locations on 
a square lattice L x L l.u. (1 < a <C L). The squares 
are fixed in space but free to overlap, and the remain- 
ing void space is filled with a fluid to which a constant, 
external force is imposed along the x axis to model the 
gravity. This system, especially at high porosities, can 
be regarded as a cross-section of a fibrous material made 
up of long, parallel fibers aligned perpendicularly to the 
flow direction. 

To make sure that the percolation threshold has its 
usual meaning, we assumed the system to be a rectangle 
of size 3L/2 x L with impenetrable walls along its longer 
side. Obstacles of size a x a were placed only in the cen- 
tral part of size L x L (see Fig. [T]). In this way any per- 
colating route through the pore space was open to flow. 
Measurements of all physical quantities were performed 
only using the data from the central, porous subsystem. 
Because our system was finite, some obstacle configura- 
tions with ip > (p c , especially close to ip c , turned out to 
block the flow completely. We rejected such configura- 
tions. We solved the flow equations numerically in the 
creeping flow regime using the Palabos (Parallel Lattice 
Boltzmann Solver) software [i§8j for a = 10 and L = 1000 
(<p c « 0.367 0). 

Figure [3] shows the tortuosity calculated from Eq. (fT2|) 
for a broad range of porosities. For 0.45 < ip < 0.9 these 
results are practically the same as those reported in [3] , 
where much smaller systems of size L x L with periodic 
boundary conditions in both directions were used and 
the values of physical parameters were extrapolated from 
those obtained for 50 < L < 300. This is consistent 
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with , where it was argued that the condition for the 
boundary and finite-size effects to be negligible in this 
model is L > 400 and a/L < 0.01. What is even more 
remarkable, very good agreement with the results of [3], 
where the tortuosity was calculated from Eq. ([9]) rather 
than from (fT2|) . indicates that the difference between the 
two formulas, resulting from how they treat reentrant 
flows, is negligible. This validates utilization of Eq. (fT2")) 
close to ip Cl where Eq. (|9]) is numerically unstable and 
hence rather useless. 

The inset in Fig. [3] depicts a log- log plot of T — 1 as a 
function of 1 — ip for large porosities (p > 0.8) at which 
the model mimics a fibrous medium. The data suggest 
that in this case T — 1 oc (1 — p) 1 with 7 = 1/2, i.e. 

T=l+ J3v /1 _ ^, (19) 

where p is a constant. This finding is at odds with most 
of conjectures about the tortuosity-porosity dependence 
for p « 1, as a vast majority of them predicts that 7 = 1. 
For example, Maxwell's formula for electrical conductiv- 
ity of a medium containing a dilute suspension of small 
spheres [3(| implies that the electrical tortuosity T c \ sat- 
isfies T c i = 1 + 5(1 - <p) as p — > 1. Similarly, Weissberg 
[3l| argued that 1 — 1 \np « 1 + |(1 — ip) is the lower 
bound for diffusional tortuosity. As for the hydraulic 
tortuosity, Mauret and Renauld in their study on fibrous 
mats assumed that T = 1 — p In ip with some constant 
p [H]. Other hypotheses include T = 1 + p(l - (p) Q, 
T = yjl +p(l - <p) [H, and T cx ip* (see 0, Hp] 
for discussion of this topic) and all imply 7 = 1. While 
the conjectures regarding diffusional or electrical tortuos- 
ity in highly porous media are well-grounded, the above- 
mentioned conjectures regarding hydraulic tortuosity are 
founded on various ad hoc approximations and even spec- 
ulations, e.g. about equivalence of the hydraulic and elec- 
trical tortuosities, and their validity is only hypothetical. 
A theoretical tortuosity-porosity relation which predicts 
7^1 was recently proposed by Ahmadi et al. [lj| . Their 
formula 

/2 P> 1 

~ Y 3 1-^(1-^)2/3 +3 

implies 7 = 2/3. The same value of 7 results from a 
formula proposed in (36|, T = p/[l — (1 — <p)%]. 

Since none of the above-mentioned formulas can be 
fitted to our numerical results, we verified that our data 
are not loaded with finite-size errors (data not shown). 
Then we investigated the flow in several highly porous 
systems. Typical examples of the velocity field in such 
systems are visualized in Fig. [1] (generated for L = 4000 
l.u., a = 10 l.u.) . As can be seen, even if obstacles oc- 
cupy only 1% of the volume so that practically each of 
them forms a separate 'island', the flow is very sinuous, as 
if restricted by some kind of solid- wall channels. These 
virtual channels are created by variations in local con- 
centration of obstacles. Since a fluid flux through a 2D 
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FIG. 4. (Color) Double logarithmic plot of the hydraulic 
tortuosity (T) as a function of a dimensionless system length 
L/a at the percolation threshold for a = 10 l.u. (symbols). 
A solid line is a fit to a power- law dependency T oc (L/a) T 
with d T = 0.19. 

channel with the no-slip boundary condition is propor- 
tional to its width cubed, the fluid passes most easily 
through the interconnected regions of low local obstacle 
concentration, whereas the regions of high local obstacle 
concentration — even if occupied by separate obstacles — 
act effectively as almost impenetrable barriers. This 
many-body effect is not present in electric, diffusional or 
inviscid fluid flows [l5| . For this reason electrical (or dif- 
fusional) tortuosity T c \ at high porosities is significantly 
lower than the hydraulic tortuosity and \dT c \/dp\ remains 
finite as ip — > 1. 

Figured] presents a log- log plot of the hydraulic tortu- 
osity dependence on a dimensionless system length L/a 
at the percolation threshold ip c w 0.367. The best fit 
to a power law yields T cx L dT with d T = 0.19 ± 0.01. 
This value is significantly exceeds the exponent d m i n = 
1.130 ± 0.002 controlling the scaling of the shortest path 
between two points on a percolating cluster [37J. This 
indicates that even at percolation most of the fluid does 
not choose the shortest-path channels. Another charac- 
teristic percolation length is the most probable traveling 
length I* , which at bond percolation scales with the sys- 
tem size as L d i with d~ t — 1.21±0.02 [38[. Using a scaling 
ansatz for the probability distribution function of a path 
length A proposed in 38] it can be shown that the average 
path length (A) ~ I* ~ L d i . Moreover, closer scrutiny of 
the method employed in [38j to generate streamlines re- 
veals that a constant-flux condition between neighboring 
streamlines was implicitly applied, just as in Eq. (|10[) . 
Hence, (A)/L oc T M , which implies 

d T = d i -l. (20) 

Our results for dT are in very good agreement with this 
conjecture. 
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V. DISCUSSION AND CONCLUSIONS 

Despite of its simplicity and common usage in vari- 
ous areas of science, the concept of tortuosity is poorly 
understood and the available literature is often mislead- 
ing, mainly because most of the theoretical research on 
this subject did not go beyond much simplified models. 
Therefore in this paper we focused on the problem of cal- 
culating the hydraulic tortuosity, defined as the average 
elongation of a streamline length in a porous medium, in 
arbitrary flow geometries. Our analysis shows that sev- 
eral existing methods of calculating hydraulic tortuosity 
differ in the interpretation of how the average stream- 
line length is to be calculated. Each of these methods, if 
applied to a system with a realistically complex geome- 
try, would yield a different tortuosity value, and only the 
method developed in [l4[ produces a number that does 
not depend on a cross-section along which measurements 
are carried out and consistently addresses the problem of 
recirculation zones. 

For incompressible fluids the method developed in [l4| 
can be reduced to calculating a ratio of the mean fluid ve- 
locity to the mean component of the fluid velocity along 
the external force direction. The two methods yield ex- 
actly the same values for regions which are connected 
by streamlines to the inlet and outlet surfaces, and dif- 
fer only in recirculation zones. As the contribution from 
recirculation zones (eddies) is expected to be negligible 
at low Reynolds number regime, both methods can be 
considered equivalent for incompressible creeping flows 
through porous media. This conclusion was confirmed by 
our numerical simulations of a 2D model of freely over- 
lapping squares. Thus, hydraulic tortuosity defined as 
the average elongation of fluid path lengths can be cal- 
culated directly from the velocity field. This not only 
greatly simplifies determination of this quantity, in ex- 
periments or numerical simulations, including complex 
3D systems, but also opens a new perspective on its phys- 
ical relevance. Many researchers doubted if an average 
path length could be defined, even conceptually, for com- 
plex flows with frequent branching and rejoining of flow 



streamlines (see [4j), but there is no doubt that the av- 
erage fluid velocity is a well-defined physical quantity. 
Moreover, the possibility of expressing the hydraulic tor- 
tuosity in terms of mean fluid velocities could be used 
to extend its definition to the case of higher Reynolds 
numbers, where the notion of individual streamlines loses 
its meaning. Note also that since no recirculation zones 
can be formed in diffusional or electrical flows, the aver- 
age elongation of streamlines in diffusional or electrical 
transport through porous media can be also exactly cal- 
culated directly from local fields, which obviates the need 
of determining individual streamlines. 

We applied the new method in two limiting cases: very 
high or very low porosities. In the former case, which cor- 
responds to a flow through fibrous materials, we found 
that the hydraulic tortuosity T scales with the porosity ip 
in accordance with T — 1 oc (1 — i^) 7 , where 7 ~ \, This 
behavior differs from that found in diffusional or electri- 
cal flows for which 7 = 1. This reflects a fact that deter- 
mination of the velocity field in a high-porosity hydro- 
dynamical system is a many-body problem, whereas the 
electric field in the same porous system can be safely ap- 
proximated as a superposition of single-obstacle solutions 
[30l l3lj . Hydraulic and diffusional (or electrical) tortu- 
osities are thus completely different quantities in highly 
porous, fibrous systems. A difference between our result 
(7 = 1/2) and a recent hypothesis by Ahmadi et al. [l9| 
(7 = 2/3) may be caused by different space dimensional- 
ity (2D vs 3D) and requires further investigations. 

When the system is at percolation, hydraulic tortu- 
osity was found to scale with the system size L as L dT 
with d T = 0.19 ± 0.01. This suggests that d T = d~ t - 1, 
where is an exponent controlling the scaling of the 
most probable traveling length at bond percolation (38j . 
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